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Multiply quantized vortices in the BCS-to-BEC evolution ol p-wave resonant Fermi gases are 
investigated theoretically. The vortex structure and the low-energy quasiparticle states are discussed, 
based on the sell-consistent calculations ol the Bogoliubov-de Gennes and gap equations. We reveal 
the direct relation between the macroscopic structure ol vortices, such as particle densities, and the 
low-lying quasiparticle state. In addition, the net angular momentum lor multiply quantized vortices 
with a vorticity k is found to be expressed by a simple equation, which reflects the chirality ol the 
Cooper pairing. Hence, the observation ol the particle density depletion and the measurement ol the 
angular momentum will provide the information on the core-bound state and p-wave superfluidity. 
Moreover, the details on the zero energy Majorana state are discussed in the vicinity ol the BCS-to- 
BEC evolution. It is demonstrated numerically that the zero energy Majorana state appears in the 
weak coupling BCS limit only when the vortex winding number is odd. There exist the k branches 
ol the core bound states for a vortex state with vorticity k, whereas only one ol them can be the 
zero energy. This zero energy state vanishes at the BCS-BEC topological phase transition, because 
ol interference between the core-bound and edge-bound states. 
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I. INTRODUCTION 

Spin-triplet p-wave supcrfluids, such as liquid 3 He, 
have brought rich physics due to interplay between the 
spin, orbital, and gauge degrees of freedom of the or- 
der parameter Various superfluid phases, so-called 
the ABM, BW, polar states and so on, arise from the 
breaking of exceptionally large symmetry group SOs(3)x 
SOi(3) x U v (l). In addition, the broken symmetry and 
boundary conditions give rise to the continuous configu- 
ration of the order parameter as non-trivial topological 
excitations, such as orbital and spin textures @, 3- It 
has also been proposed [4, 5] that p-wave superfluids can 
exhibit another type of the phase transition, which does 
not involve symmetry breaking but changes a topological 
charge. The topological phase transition is driven by the 
enhancement of the effective pair interaction. 

Recently, p-wave Fcshbach resonances on 6 Li and 40 K 
atoms have been observed by sweeping magnetic field in 
experiments The Fcshbach resonance of colliding 

atoms into an £-wave bound state allows one to manipu- 
late the inter-atomic interaction in ^-wave channel, from 
the regime of weakly interacting atoms to tightly bind- 
ing molecules regime. In the weakly interacting regime, 
the Fermi system forms Cooper pairs, where the low en- 
ergy quasiparticle excitation is characterized by a p-wave 
Cooper pair potential. In contrast, the p-wavc molecules 
undergo Bosc-Einstein condensation (BEC), which in- 
volves the isotropic excitation gap uniquely determined 
by the binding energy of molecules. It has been clari- 
fied that in contrast to the BCS-BEC crossover which 



appears in s-wave supcrfluids [111 ], these two regimes are 
not smoothly connected but give rise to the topological 
phase transition 0, [H, EH, EH • The experimental obser- 
vation of p-wave Feshbach resonances is the first step to- 
wards the further realization of the p-wave BCS-to-BEC 
topological phase transition. 

Here, it is natural to expect that the extremely strong 
coupling regime of p-wave superfluids contains intrigu- 
ing vortex structure, because the coherence length of 
the order parameter £ becomes comparable to the mean 
interparticle distance fcp 1 . The phase of the order pa- 
rameter rotates by 2t:k around a quantum vortex, where 
KgZ because of its single valuedncss. Hence, the quasi- 
particles traveling across the vortex core experience the 
abrupt shift of the phase, leading to the low-energy An- 
dreev bound state [14[ or the so-called Caroli-de Gennes- 
Matricon (CdGM) state [H|. In the regime with fc F 1 
of s- wave sup erfluids, it has been first demonstrated in 
Refs. [H, U3 that the energy levels of the CdGM state 
become discrete, leading to the strong depletion of the 
particle density around the core. The analysis on the 
quantum depletion has been theoretically extended to the 
BCS-BEC crossover regime (l8l - l23| and the depletions in 
the particle density were experimentally observed in ro- 
tating Fermi gases with an s-wave resonance [241 . [25j . 
Furthermore, the structures of giant vortices with |re| >1 
have been extensively studied in s-wave supcrfluids by a 
large number of authors [26H32| . 
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Turning now to p-wavc supcrfluids, the CdGM states 
are found to be intriguing in the following two senses: (i) 
The low energy states in the weak coupling BCS phase 
are topologically distinct from them in the strong cou- 
pling BEC phase @, Q3]. (ii) The lowest eigenenergy 
may be exactly zero in the BCS limit that the chemi- 
cal potential fi is equal to the Fermi energy E~p [33 38J. 
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The zero energy state (ZES) bound at the vortex core 
and at the edge can be described as the Jackiw-Rebbi 
solution in the one-dimensional Dirac equation at the do- 
main wall }39l - l44| . The remarkable feature of the ZES is 
linked to the fact that the creation operator is identical 
to its own annihilation, called the Majorana fermion [46|. 
In addition, since the host vortices associated with the 
ZES obey neither Fermi nor Bose statistics, called the 
non-abelian statistics [H, H?J , this system will offer the 
promising method of the fault-tolerant quantum compu- 
tation |48l - [50l | . It is known that they can be observed in 
spin-polarized p-wavc superfluids with the time-reversal 
symmetry breaking, e.g., the k x ±ik y chiral state. Hence, 
p-wave resonant Fermi gases can be base platform to 
check various theoretical issues because of the high con- 
trollability, such as two dimensionality [l2|, HH • 

The aim of this paper is to clarify the macroscopic 
structures of giant vortices from the fully microscopic 
point of view in the BCS-BEC evolution of p-wave res- 
onant Fermi gases. In the previous paper [52|, we re- 
vealed the roles of the ZES for visualization of p-wave 
superfluidity through the particle density depletion and 
net angular momentum. In this work, we expand this ar- 
gument into other vortices with winding numbers \k\ > 1. 
Here, we use the fully microscopic theory based on the 
Bogoliubov-de Gennes (BdG) equation, which enables 
the systematic study in the BCS-to-BEC regime. It is 
demonstrated that the low-energy quasiparticle spectra 
is reflected through the quantum depletion of the parti- 
cle density at the core. In addition, we calculate the net 
angular momentum, which is found to provide direct ev- 
idence for p-wave superfluidity. Throughout this paper, 
we focus on the zero temperature limit because for finite 
temperature regimes the pairing fluctuation effect which 
is not taken into account here becomes important (53j . 

Furthermore, we discuss the existence of the zero 
energy Majorana states in the BCS-to-BEC evolution 
regime. It has been revealed in the BCS limit [43| and 
the more generic situation (45j of p-wave superfluids that 
there is only a single zero energy Majorana state for odd 
vorticity and none for even vorticity. This is contrast to 
the index theorem for zero energy eigenstates of the rela- 
tivistic Dirac Hamiltonian [54|, [55| and the quasiclassical 
analysis of the p-wave BdG equation (3, |35| ■ Here, we 
reproduce the prediction and further show that the ZES 
in an odd-vorticity vortex disappears in the BCS-BEC 
topological phase transition because of the quasiparticle 
tunneling between the core- and edge-bound ZES's. It 
is also confirmed numerically that their wave functions 
are characterized by the modified Bessel function in the 
vicinity of the topological phase transition [45[ . 

In the following section, we describe the details on 
the self-consistent calculation based on the BdG and gap 
equations. In this section, we summarize the results of 
the vortex-free state in the BCS-to-BEC evolution. The 
vortex structures, such as the profiles of the order pa- 
rameter and the particle density, are presented in Sec. III. 
These results are discussed with the knowledge of the low- 



energy quasiparticle structure. Furthermore, in Sec. IV, 
we demonstrate that the lowest eigenenergy becomes zero 
when the winding number is odd, whereas it lifts from 
zero in the vicinity of the BCS-BEC evolution. The fi- 
nal section is devoted to conclusion and discussion. The 
detailed derivations of the self-consistent equations and 
the CdGM states with arbitrary winding number are in- 
cluded in Appendices. 



II. THEORETICAL FORMULATION 

A. Self-consistent equations 

Here, we start with the mean-field Hamiltonian of spin- 
less fermions with the mass M and the Nambu spinor 
*(ri) = [V(ri),V t ('l)] T , 

H = E + ^Jdr 1 j dra^nj^n.ra^ra), (I) 

where ijy and tp are the creation and annihilation opera- 
tors of fermions. The matrix JC is given as 



with 



H (r)S(rx -r 2 ) A(ri,r 2 ) 
A*(r 2 ,n) -HS(n)S(n -r 2 ) 



Ho(r) 



2M 



- jji + ift(xd y - yd x 



(3) 



Throughout this paper, we set h=kg = 1. The rotation 
frequency Q is set to be zero throughout this paper. The 
pair potential A(ri,r 2 ) is defined as 



A(n,r 2 ) = -~V(r 1 ,r 2 )(ip(ri)^(r2)}. 



(4) 



Under p-wave pair potentials A m (r) (m = 0, ±1), 
the quasiparticle eigenstate with the wave function 
[u v (r),v u {r)] T is described by the BdG equation 



H (r) n(r) 
-H*(r) -HS(r) 

U m=0,±l 



u v {r) 
v v (r) 



A m (r)7>„ 



EL 



v v {r) 



-V m A m (r) 



, (5a) 



, (5b) 



where we introduce the p-wave operators V±i = 
=F (d x ± id y ) and Vo = d z . The details on the deriva- 
tion are described in Appendix A. The wave functions 
must satisfy the orthonormal condition 



[ut(r) Ufl (r) + u*(r)u M (r)] dr = 5„ tl 



(6) 



The pair potential is self-consistently determined by 



(7a) 
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Ao(r) = - 



(7b) 



where $(ri,r2) is defined in Eq. (|A9b[) with the eigen- 
states of Eq. ([5]). 

The BdG equation ([5]) and gap equation ([7|) provide 
a qualitative formalism for the BCS-to-BEC evolution 
regime Il2l . [ill ] , by combining with the following num- 
ber equation and by using the renormalized coupling con- 
stant. First, the total particle number N is fixed by ad- 
justing the chemical potential /i, where the total number 
is expressed as N~ J p(r)dr with the particle density, 



(8) 



The sum denotes the summation for all the eigen- 
states with the positive and negative eigenvalues. f(E) = 
l/(e B / T +l) is the Fermi distribution function at temper- 
ature T. In this paper, we set temperature to be T = 0. 

The gap equation ([7]) involves two divergence terms 
proportional to E c and In E c [l2|, [5(| [57} where E c is the 
cutoff energy. The dominant divergence can be removed 
by replacing the bare coupling constant to the renormal- 
ized one 



|iW 



s fc 2efc 



(9) 



with the volume of the system S, eu = k 2 /2M, and 
T(k) = k/k - Here, E\> is an eigenvalue of the Schrbdinger 
equation for two fermions interacting via the pairing po- 
tential V [56|, H3] ■ Eh is real and regarded as the two- 
body bound state energy in vacuum when E^ is negative, 
while it has an imaginary part for positive Ey,. However, 
it is known that this imaginary part is negligible in the 
vicinity of a p-wave resonance |12j . Hence, the coupling 
constant g m parametrized by E^ remains real and neg- 
ative for all values of Eb, which can remove the leading 
term of the ultraviolet divergence in gap equation (JT]). 
Note that in the case of a two-dimensional geometry, the 
resulting gap equation ([7]) with Eq. © still contains the 
logarithmic divergence on E c . 



B. Cylindrically symmetric system 

In order to study vortex structures, it is convenient 
to introduce the cylindrical coordinate r = (r, z). The 
k x ± iky pairs have the phase winding ±1, since V± is 
expressed as V±\ = ^fe ±ie (d r ± ^dg). We assume the 
cylindrical symmetry of the order parameter 

A m (r) = A m (r)e^ e , (10) 

which describes the quantized vortex with a winding 
number K m 6Z centered at the origin r = 0. 



The quasiparticle wave function is then expressed in 
the cylindrical coordinate as 



u u (r) 




_ u v (r) _ 





(11) 



Here, we impose the periodic boundary condition with 
the period Z along the z-axis, i.e. , q = 2im z / Z with n,eZ 
and the axial symmetry leads to £ u ,£v& Z. Since we here 
consider a two-dimensional system, however, the axial 
quantum number along the z-axis, q, is set to be q = 0. 
From the BdG equation we find the condition that u v 
couples with v v through the following angular momentum 
relation, 



£ v + k +1 + 1 = £, 



(12) 



with the azimuthal quantum number £ € Z. Also, one 
obtains the condition for the phase winding of A m (r) as 



K + l — — 1 — K— \ 



(13) 



In addition, the following symmetry relation for eigen- 
states with q = or with Ao = should be noted: 
The positive energy state Eg and [m£,^] t with £ is 
symmetric to the negative energy state — E-£ +K+ i and 
[ue,v e ] = [v*_ e+K+1 ,u*_ e+K+1 ] T with -£ + k+1. 



C. Calculated system 

In the current work, we focus on a two-dimensional 
system where the Yi^(k) orbital state of the pair po- 
tential is neglected, i.e., Ao(r) = 0. The rigid wall 
boundary condition is imposed on the wave functions, 
u v (r = R) = v u (r = R) = 0, at R = 50k p 1 . In this pa- 
per, we consider the strong coupling regime nearby the 
resonance, in which the coherence length, 



fep 



MA 



Ml, 



(14) 



is comparable to the average of the interparticle spacing, 
i.e., £ ~ kp 1 . Hence, the radius R which we set here 
is enough to recover the pairing field from the vortex 
center, i?^>£. The resulting system consists of the A^ = 
600 fermions in a single hyperfine spin state. All the 
quantities are scaled by using the length unit fcp 1 and 
the energy unit E-p, where Ep = k F /2M is the Fermi 
energy. The energy cutoff is set to be E c = 6QEp. We 
also assume ko = k-p. 

In the vortex-free state of two-dimensional spinless 
Fermi systems, the k x ± ik y pairing states are energeti- 
cally degenerate, called the orbital ferromagnetic or chi- 
ral state. The non-zero vorticity, however, breaks the de- 
generacy and vortex states with a positive winding num- 
ber becomes distinguishable from the negative winding 
states. Hereafter, without the loss of generality, we con- 
sider the system that A+i is dominant. 
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FIG. 1: (Color online) (a) Chemical potential /iasa function 
of i?b in the vortex- free state with re = 0. The maximum value 
of the pair potential A/E-p — max \A+i(r)\/Ep is plotted as 
functions of £ [—3, 4] in the inset of (a), (b) Bulk excitation 
gap in the vortex- free state as a function of Eh- 



D. Vortex- free state 

First of all, we present in Fig. Q] the numerical re- 
sults of the self-consistent equations ([5]) and ([7]) with 
Eq. © in the vortex-free state (re = 0). The chemi- 
cal potential fj, and the maximum value of the pair po- 
tential A = max|A+i(r)| arc plotted as a function of 
E\, in Fig. HJa). Feshbach resonance allows one to ma- 
nipulate the effective coupling constant g±i, where the 
Cooper pairs realized within (i > 0, say the BCS phase, 
turn to molecular bosons in the BEC phase when fx < 
@; EH EH • The coupling constant g m is parameterized 
with E\j through Eq. In E\> > 0, \i approaches the 
Fermi energy Ep with increasing E\>. It is seen from the 
inset of Fig. [3a) that throughout the whole range of Eb, 
the amplitude of the pair potential A yields no singular 
behavior. 

Figure [Hb) shows the low-energy quasiparticle spec- 
trum in the bulk. In the BCS limit E^^Ep, the lowest 
excitation gap in the bulk is characterized by the disso- 
ciation energy of the Cooper pair min \E\ = A. In the 
vicinity of \fj.\ <C Ep, however, the energy gap turns to 
min l-El = \fj.\, which indicates that the excitation becomes 
gaplcss. In our calculation, this occurs at E^/Ep sa — 1.0, 
as seen in Fig. [TJb) . Read and Green [j| predicted that 
the point with /i = separates two topologically distin- 
guishable phases, such as the BCS (Eh/Ep > — 1.0) and 
BEC (Eb/Ep < — 1.0) phases. This is contrast to the case 
of the s-wave pairing; The excitation spectrum in the 
fi > regime has a gapful with |A|, which continuously 
turns to min \E\ = -y/|/i| 2 + |A| 2 as fi becomes negative 
across fi = 0, i.e., the BCS-to-BEC crossover |11|. 



III. VORTEX STRUCTURE 
A. Pair potential 

In two-dimensional systems, the order parameter re- 
duces to two components, A±i, because of q = 0. Here, 
we consider vortex states that the majority component 
A_|_i has a winding number |re| < 3. As seen in Eq. (|13l) . 
the axial symmetry of the system requires the winding 
number to satisfy the condition re+i = k—i — 2 = n. Hence, 
we take account of 6 combinations of winding numbers, 
(k+i,k-i} = <-3,-1>> (-2,0), (-1,1), (1,3), (2,4), and 
(3, 5). Here, we shall discuss all the vortex states, except 
for (k+i, K-i) = (—1,1), because its winding configura- 
tion was studied in our previous work [52[ . 

Figure [5] shows the profile of the pair potentials A±i(r) 
near the vortex center r = in the k — ~ 2, —3, 1, 2, 3 vor- 
tex state. Here, it is convenient to categorize the vor- 
tex states in terms of the sign of K. The "positive" 
("negative") vortex state is defined as the state that 
the vorticity re is parallel (anti-parallel) to the chiral- 
ity of the orbital motion k x + ik y of pairs. In general, 
a quantum vortex is defined as the phase singularity at 
which the amplitude of the order parameter vanishes as 
A(r = 0) = 0. The order parameter recovers to the con- 
stant value within the coherence length £. The slope of 
A(r) near r = is characterized by the winding number, 
i.e., lim r _>o A m (r) oc r' Km l . It is seen in Fig. [21a) that in 



(a) u.6 




FIG. 2: (a) Profiles of the pair potential A+i(r) in re = +1 
(solid line), +2 (dashed line), +3 (dotted line) at E h /E F = 1.0. 
(b,c) Profiles of pair potentials A+i(r) (solid line) and A_i(r) 
(dashed line): (b) re = -2 and (c) re = -3 at E h /E F = 1.0. The 
inset in (a) shows the vortex radius £ C orc normalized by fcj 
as a function of Eh- The solid and dashed lines denote f CO re 
and £/2, respectively. 
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the case of the positive vortex state, the winding number 
K m determines the size of the vortex core, which becomes 
larger as K m increases. The vortex core of A+i in the case 
of the positive vortex state is always smaller than that 
in the induced component A_i, because of k+i < 
Hence, the giant cores of A+i leave empty and the result- 
ing pair potential is effectively describable with the single 
component A + i(r), where A_i(r) with the larger vortex 
core is suppressed. In addition, the core size of the n = +l 
vortex can be quantified as £^ e = linv-)0 |A +1 (r)|/r with 
A +1 (Y) = A +1 (r)/A [l6|. As seen in the inset of Fig. [^a), 
the strong coupling effect gradually shrinks the core ra- 
dius £ C orc, where the vortex core size is of the order of 
the atomic scale ps k^ . 

In contrast, since |k+i| > in the case of the nega- 
tive vortex states with k + i < — 1, the induced component 
A_i always has the vortex core smaller than that of A+i. 
It is seen in Figs.[3Jb) and^c) that the A_i component 
is induced inside of the vortex core of majority A+i(r). 
In particular, the vortex core of A+i with k= —2 is filled 
in by the induced component A_i which is the vortex- 
free state with K-i = 0. Hence, the quasiparticlc struc- 
ture at the core exhibits gapful. In the same sense, since 
the vortex structure of k= — 3 in the vicinity of the origin 
is dominated by the induced component, the quasiparti- 
clc structure and particle density at the core arc found 
to be almost same as those in the single- vortex state of 
A_i with K-i= — 1. 

Note that these vortex structures and the local den- 
sity of states are also discussed in Ref. [HH within the 
quasiclassical approximation, corresponding to the BCS 
limit with fi = Ep. In the BCS regime, e.g., E^/E-p = 1.0, 
the pair potentials of all the vortex configurations which 
we obtain here is found to be consistent with those in 
Ref. HI. 



B. Low- lying spectra 

We draw attention to the quasiparticlc excitation spec- 
trum under the pair potential displayed in Fig. [21 In 
particular, we focus on the vortex states with positive 
winding numbers k = +1, +2, +3. In Fig. [3l we plot 
the low-lying energy spectra of their vortex states at 
three regimes: The weak coupling BCS regime {E\ } /E-p = 
1.0), the vicinity of the BCS-to-BEC transition point 
(E h /E F = -0.6), and the BEC regime (E h /E F = -1.2). It 
is found that two branches are embedded inside the bulk 
energy gap \E\ <0.4Ep in the BCS regime of k=1. The 
low- lying branch labeled as the "Edge" in Fig. [3] consists 
of the eigenstates whose wave functions are localized in 
the edge of the system. This edge mode branch exists in 
the other vortex state, e.g., n = 2 and 3, as seen in Fig. [3] 
The eigenenergy of this edge mode in the k x + ik y pairing 
state with a vortex winding number k obeys the linear 
dispersion on the azimuthal quantum number £ |44| . 



E f = - £- 



1 



where the energy spacing is found to be e w A/kpR = 
O.OOSE'p at E h /E F = 1.0. Note that the eigenenergy of 
the edge mode can become zero at I = , when k is 
odd. 

In the weak coupling regime (E^/Ep = 1.0), the eigen- 
state belonging to the other branch inside the bulk excita- 
tion gap -E = ±Ass0.4£f is identified as the core bound 
state, i.e., the CdGM state labeled as the "CdGM" in 
Fig. [3] The eigenenergy of the CdGM branch in the 
k x + iky pairing system with arbitrary winding number 
Ky^O can be analytically expressed with n £ Z and q = 
as 



El t n — — 



1 



w 



1 



(16) 



This is obtained from the BdG equation ([5]) by extending 
the procedure by Caroli et al. [15| to the chiral p-wave 
system, which is valid for the eigenstate within \£\ <C&f£- 
The expressions of two coefficients wo,i are described in 
Eq. (|B19j) . As seen in Eqs. and (|B19|) . even though 
the minority component A_i(r) contributes to the coef- 
ficients LJo.i, the eigenenergy in Eq. (|16p maintains the 
zero energy state with £ = ^4^- when k is odd. Since 

these are approximated as an( l wi«|A in the 

BCS limit (A<CEp), one reads wo<wi- The numerical 
results in Fig. [3] reproduce the linear behavior on £. The 
wave function is exponentially localized inside the core: 



ut(r) 
vt{r) 



fe,n(k+ r ) 
//-„-i,n(fc-r)e-*C"+i)« 



(17) 



where M is the normalization constant and 



ft,n(k±r) J,(fc±r)e-^ / o" [A+l(r ' ) - A - l(r ' )]rfr '. (18) 



kfj, for the 



(15) 



Here, k± = k^± MEi, n /k^. Note that k± 
low-lying eigenstates in the BCS regime. The further de- 
tails are described in Appendix B. As the strong coupling 
regime is approached, however, the CdGM wave function 
spreads over the whole region of the system. This will be 
discussed in Sec. IV. 

It is found that the number of the CdGM branches 
uniquely depends on the winding number k in the weak 
coupling BCS regime, e.g., three CdGM branches ap- 
pear in the k — +3 vortex state. In addition, the low- 
est eigenenergy of both the CdGM and edge states in 
Eqs. (|15p and (fT6"| can be zero only when k is odd, which 
is numerically demonstrated in Fi g. Pp This reproduces 
the prediction derived in Refs. (43l |45|| that there is only 
single zero energy state for odd vorticity and none for 
even k. It is obvious that the CdGM wave function with 
£ exhibits ug(r) oc Je{k^r) m and Vt{r) ~ r \i-K—i\ 
at r — > 0. Hence, the asymptotic wave functions of 
the zero energy states with £ = (k + l)/2 E Z exhibits 
ug(r) =vp(r) «r' K+1 l/ 2 . The further details on the zero 
energy states will be described in Sec. IV. 

In the vicinity of the BCS-BEC phase transition <C 
E F , e.g., E\j/E-p = —0.6, the bulk energy gap is char- 
acterized as min |i?buik| = as seen in Fig. [TJ In this 
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FIG. 3: (Color online) Quasiparticle excitation spectra at E^/Ep — 1.0 (left column), —0.6 (center column), and —1.2 (right 
column) in the vortex states with k = 1 (top row), k = 2 (middle row), and k = 3 (bottom row). "CdGM", "Edge", and "ZES" 
denote the branches of the CdGM, edge, and zero energy state, respectively. 



regime, the amplitude of the pair potential is compara- 
ble to the Fermi energy, A sa Ep . This means that the 
energy level distance of the CdGM states becomes com- 
parable to the Fermi energy, i.e., uiq ps Ep in Eq. (fT6"|) . 
In contrast, the bulk excitation gap, i.e., the lowest en- 
ergy of the continuum states, becomes proportional to 
the amplitude of the chemical potential, min |i?buik| = 
Since \fj,\ < Ep in this regime, the CdGM states are 
merged to the continuum state with the bulk excitation 
gap min |-Ebuik| = < uj w A 2 /Ep. Hence, as seen in 
Eb/Ep = —0.6 of Fig. [3l the low-lying spectra in this 
regime consist of the branch of the edge state with or 
without the zero energy state and bulk excitation gap. 
Note that the bulk excitation gap is still larger than 
the level distance of the edge mode e in Eq. (|15j) which 
is determined by the inverse of the radius of the system. 

Beyond the topological phase transition E^/Ep = — 1.0, 
the low-lying spectra becomes trivial, where the excita- 
tion gap is uniquely characterized by min \E n \ = It 
is seen from Fig. [3] that the branch of the edge state and 
the zero energy state disappear in the BEC phase. 

However, it should be emphasized that even in the 
BEC phase, the spectra is asymmetric with respect to 
I. This arises from the fact that the BdG equation ([5]) 
still requires the one-to-one mapping between [u£,ve] T 



l+K+ii v -e+K+i 



] with — _E_£ +K _|_i. This 



with Ei and [i 

asymmetry on I gives rise to the non-trivial angular mo- 
mentum. The total angular momentum per particle is 



estimated by the fully self-consistent calculation as 

(L z ) _K + 1 
HN ~ 2 ' 



(19) 



in the whole range of E\> in the k x + ik y pairing state, as 
seen in Fig. [4] The deviation of the self-consistently ob- 
tained (L z ) from Eq. (|19p arises from the vortex winding 
and chirality of the induced component A_i(r). In the 
weak coupling BCS limit, the total angular momentum 
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FIG. 4: (Color online) Total angular momentum (L z ) 
as a function of E^ in various vortex states with n = 
-1,0, +1, +2, +3. The dashed lines denote Eq. 1(15)1. 
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per a particle (L z )/N consists of the angular momentum 
due to the vorticity hn/2 in addition to the chirality of the 
k x + iky channel +H/2, where the latter originates from 
the linear dispersion of the edge state in the BCS regime 
plil [59j . Figure [4] demonstrates that this argument can 
be expanded to giant vortex states in the vicinity of the 
BCS-BEC evolution. Note that since (L z )/N &Hk/2 in 
s-wave supcrfluids, the measurement of the total angular 
momentum through the shift of quadrupole frequencies 
provides direct evidence for p-wave superfluidity in ex- 
periments [52| . 



C. Particle density depletion and local density of 
states at vortex cores 

The direct observation of the low-lying quasiparticlc 
spectra is experimentally challenging problem. Neverthe- 
less, it has been discussed in s-wave superfluids with arbi- 
trary winding number (l7 L 18 L 20l-|23i l28l - l32j and p- wave 
superfluids with |«| < 1 52l. 60| that the local particle 
density p(r) around the vortex core reflects the low-lying 
CdGM spectrum. 

We here extend this analysis to giant vortex states in 
p-wave BCS-BEC evolution. Among the various possible 
vortices, the particle density depletion in the ft = ±l vor- 
tex state was discussed in the weak coupling regime [6(J ■ 
In addition, the re = — 1 vortex in the vicinity of the BCS- 
to-BEC evolution was studied in Ref. [EH- It is found 
that the core of the k = — 1 vortex in the weak coupling 
BCS regime 3> Ep is invisible through the density pro- 
file, which reflects the fact that the vortex core is filled in 
by the zero energy CdGM state with £ = 0, i.e., the Majo- 
rana state. As the BCS-BEC topological phase transition 
point is approached, however, since the zero energy state 
lifts from zero to finite energy, the vortex gradually be- 
comes visible. As a result, the quantum depiction of the 
k = — 1 vortex turns out to be closely associated with the 



instability of the zero energy state with £ = across the 
BCS-BEC topological phase transition point. 

Let us now attempt to expand this argument into the 
vortex states with the other n. Figure [5] shows the par- 
ticle density profile around the giant vortex cores with 
k = +1, +2, and +3. In the weak coupling limit E^^> Ep 
of all the vortex states, the vortex core region is filled in 
by particles, while p(r = 0) becomes zero in the BEC limit 
Eb/Ep <^—l. In the k = +1, the particle density at the 
core p(r = 0)/ max p(r) is plotted in the inset of Fig. [5ja) , 
which shows smooth suppression toward zero as E\, de- 
creases. Here, the core of k = +1 vortex is quantitatively 
more visible than that of the k = — 1 vortex in the whole 
range of E^. In the giant vortex states with k = +2 and 
+3, however, it is seen from the inset of Figs. [SJb) and 
[SJc) that discontinuity appears in the depletion at cer- 
tain E h , e.g., E h /Ep^0.5 for k = +2 and E h /Epm-0.1 
for k = +3. 

These behaviors on the particle density, such as the 
depletion and its discontinuity, are closely linked to the 



CdGM state with £ = 0. In the BCS regime, we have 
seen in the previous section that the low-energy exci- 
tation spectra consist of the CdGM and edge states in 
addition to the bulk excitation gap. Hence, it is nat- 
ural to decompose the particle density around the core 
into two contributions: p(r) = Pc6Gm(t) + pbuik( r )- Here 
Pbuik(r) = E|s|>a \ u d r )\ 2 f( E ) ma y bc assumed to be 
uniform in this regime, while one finds 



PcdGM(r^0)= J2 M^0)| 2 « 



(20) 



<>«±1 
c — 2 



with Eqs. (TT6) and (flT)). One finds that pcdGM(r-^O): 



E 



00 1 



"° for the K=— 1 vortex and /?cdGM(( r — >-0)« 
£4 / ■ ■ ~ r 2 for the k = +1. Hence, this implies that 
the quantum depletion displayed in the inset of Fig. [5ta) 
reflects the discrepancy of pcdGM<V = 0) between k = ±1 . 

Following the same argument, the giant vortices with 
the larger winding number are also understandable. In 
the weak coupling limit Aq/Ep -C 1, it is found from 
Eq. (TO)]) that two (three) CdGM states with £ = can 
exist inside the bulk gap in the k = +2 (+3) vortex state. 
For instance, E^q rj ±-j + C(A 2 / Ep) is obtained from 
Eq. (fT6|) in the case of k = +2. Only the negative energy 
state contributes to pcdGM( r = 0) in Eq. (f2"0"T) . However, 
it is predicted that the correction with C(A 2 /Ep) shifts 
the energy E^q upward as the strong coupling regime 
is approached. As a result, the shift of E^ from the 
negative to positive energy leads to the abrupt depletion 
of p(rwO). This argument is also applicable to the case 




FIG. 5: (Color online) Spatial profiles of the particle den- 
sity p(r) in k = +1 vortex state (a), n — +2 (b), and 
« = +3 (c). The insets in (a)-(c) show the density deple- 
tion p(r = 0)/ max p(r) at the vortex core as a function of Eh- 
For comparison, we plot p(r = Q) in the k = — 1 vortex state 
with open triangles in the inset of (a). 
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of k = 3, where the three lowest eigenenergics of the £ = 
CdGM state are found to be E^q as 0(A 2 /E f ) and 
E^o as ±A + 0(A 2 /^ F ). Only the distinction between 
k = +2 and +3 arises from the quantitative difference 
of Eg=o, resulting in the fact that the critical value of 
Eh, at which p(r — 0) exhibits jump seen in the inset of 
Figs.[SJb) andJSJc), depends on k. 

To further clarify the relation between p(r = 0) and 
the low-lying quasiparticle state, we present in Fig. [6] 
the local density of states (LDOS) at the vortex core 
N{r = 0,E) and at the bulk Af(r = §,E) in various 
vortex states, where the LDOS is defined as 



(21) 



with the Lorentzian function 5 v (z) = (r//2) 2 /[z 2 + (r]/2) 2 ] 
(see for instance, Ref. [6l|). In Fig. [6j we set 7y = 0.02i?F- 
In the case of k= +1, as we have expected, it is seen 
in Fig. EJa) that Af(r = 0, E) has the sharp peak about 



R/2, E) 




E/E s 



FIG. 6: (Color online) Local density of states at r — (left 
column) and r = R/2 = 25/cf (right column) in the range of 
Eb/Ep G [—1.0, 1.0]: (a) k = +1, (b) k = +2, and (c) k = +3. 
The arrows in the panels denote the CdGM states with £ = 0. 
For the LDOS in the case of k=—1, see Ref. 1521. 



E/E-p as 0.2 at Eh/E-p = 1.0 denoted by the arrow in 
Fig.jBfa). This peak reflects the lowest CdGM state with 
£ = 0, which is not occupied and leads to pcdGM (r = 0) = 0. 
As the strong coupling regime is approached, the in- 
tensity of the peak around E/E-p as 0.2 is suppressed 
and merges to the other bulk excitations. In addition, 
Pbuik(?~ = 0) is constructed from the eigenstates around 
E/Ep as — 0.5. The negative energy states gradually be- 
comes empty in the vicinity of the BCS-BEC evolution 
(E h /E F = -1.0), leading to p(r = 0) = 0. In the other 
vortex states with n = +2 and +3, the peak originating 
from the £ = CdGM state shifts to the positive energy 
region at E^/Ep as 0.5 and —0.1, respectively. This re- 
sults in the abrupt jump of the particle density depletion 
p(r = 0) in the inset of Fig. [5] 



IV. ZERO ENERGY MAJORANA STATES 



It is found from Eqs. ([15]) and (|T6| that the energy of 
the lowest CdGM and edge states can be zero when k is 
odd and the weak coupling limit A<g^Ep is approached. 
In contrast, the low-energy spectrum yields an isotropic 
gap uniquely determined by \p\ in the BEC phase with 
M<0. 

Figure [7] shows the lowest eigenenergies in the vicin- 
ity of the BCS-BEC evolution E h /E F G [-1.8,0.0]. In 
the BEC regime E-^/Ep = —1.8, all the vortex states 
have the energy gap comparable to \p\ =0.21£'f, regard- 
less of k. In the case of k = 2, the lowest eigenener- 
gies stay around E/Ep «0.01 even if the weak coupling 
BCS regime is approached beyond the transition point 
E\j I Ep = — 1.0. These turn out to be the lowest edge state 
with £ = +l, whose energy is in good agreement with the 
analytic result in Eq. {TSJ}, E^ +l = A/(2k F R) as 0.01£ F 
with i? = 50fc F " 1 and AasEp. 




-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 
E h /E» 

FIG. 7: (Color online) Lowest eigenenergies with logarithmic 
scale as a function of Eb/EF G [—1.8, 0] in vortex states with 
K = — 1 (open circles), +1 (filled circles), +2 (triangles), and 
+3 (squares). The solid line corresponds to \p\. 
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In contrast to the even k case, the eigenenergies expo- 
nentially shift to zero when k is odd. These states are in 
a consequence of the quasiparticlc tunneling between the 
edge and CdGM states with £=(k+ 1)/2<gZ. Assuming 
the limit of R— >oo, it is known that each state can have 
the exactly zero energy, as seen in Eqs. (fT5|) and ([To]) with 
1=(k+ 1)/2gZ and Refs 0, EH . It is known that the 
zero energy eigenstates exhibit the Majorana property, 
such as r t _ 1 = r £ y ) , resulting from UE=o(r) = v^ =0 (r) 

As a consequence of finitencss of the system with 
kpR = 50, however, these two Majorana wave func- 
tions bound at the core u c t {r) and at the edge w|(t") 
are hybridized with each other, resulting in the sym- 
metric and anti-symmetric states, u\ = [u'j, + m^]/\/2 and 
u e = i u e ~ u l]/v2- The hybridization leads to the split- 
ting of the degenerate zero energy states to the finite 
energy ±E. Based on the argument about the symmetry 
of the eigenstates in Sec. II B, if the wave function of the 
state with +Ei is [ui,vi] = [ug,uf\, then one finds that 
the other is [uf*,uf] for —Eg. 

Figure (5Ja) shows the wave function |ufci(r)| and 
|v£=i(V)| of the lowest energy state with Eg^/Ep = 
— 1.9 x 10 -5 in the BCS regime of the k = +1 vortex 
state. Here, the eigenfunction of this state turns out to 
be = [u\,uf\. It is seen that the wave function 

is exponentially localized around the core {r = 0) and 
around the edge (r = 50fc F ~ 1 ) within the coherence length 
C = 2fc F /(A/A)=4.8fc F 1 with A = 0.42^. Note that due 
to I 7^ 0, |uf(r)| = |i>f(r)| w r is obtained at ) — > 0. The 
period of rapid oscillation is found to be an order of the 
27r/fc M = 7.2fcp 1 with fi = Q.757Ep, resulting from Jg{k^r) 
in Eq. flT7|). 

These results displayed in Figs. [7] and [8] arc consistent 
with the analytic expression in two-vortex systems (62l | 
and the numerical results in vortex-antivortex systems 
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FIG. 8: Wave functions \ue(r)\ (solid curve) and |i>£(r)| 
(dashed curve) of the CdGM state with £ = at E b /E F = 2.0 
(a) and —1.0 (b). The eigenenergies are E/Ef = — 1.9 x 10~ 5 
in (a) and E/E F = -0.024 in (b). 



[63, |64[ and in plural- vortex systems 65 1. In addition, 
the zero energy states are found to exist in chiral p-wave 
systems within tight binding approximation [66|, • It 
was proposed in Ref. [62| that due to the rapid oscillation 
in the wave functions, the Friedel-like oscillation of the 
splitting energies appears. This is numerically demon- 
strated in Ref. [65(. Note that the splitting of the zero 
energy Majorana states, which is critical to the deco- 
herence in the topological quantum computation, is also 
observed in other systems, such as the non-abelian quasi- 
holes of the v = | fractional quantum Hall state (68l [1^) , 
Kitaev's honeycomb lattice model [70l| . and the generic 



anyon model 71 1 



As E\, decreases from Ef to the negative value, the 
chemical potential touches to the zero, e.g., ^l^O.OIEf 
at Eb/EF = —1.0. Then, the length scale of the wave 
functions ut and vg spreads over the whole system, e.g., 
2n/k tl = 65k F 1 >Ra,t E h /E F = -1.0, as seen in Fig.^b). 
This also reflects the topological phase transition at 
Eb/EF ~ —1.0 with /i = 0, where the bulk excitation 
becomes gapless @. 

Finally, in Fig. [9J we compare our numerical results 
for the wave functions \ug\ of the lowest CdGM branch 
of K=+l and +3 vortices with the analytic expressions 
described in Eq. (JT7J) and in Refs. [HEI- As described 
in Appendix B, the CdGM solution assumes the weak 
coupling limit /cp£3>l, i.e., A<C-Ep, resulting from the 
linearization of a fcrmion dispersion on the Fermi sur- 
face. The spatial shape of the resulting wave function 
consists of two different length scales, such as J^(fc M r) and 



The former describes 
few 1 and the latter is 



e fc F A( - r ^ dr . as seen in Eq. (fTSI) . 
the shorter length scale with fc~ 1 R 
slow function varying over ^^fcp 1 . The alternative solu- 
tion for the zero energy solution of the k — — 1 vortex was 
discussed by Tcwari et al. [z3 in both the BCS and BEC 
regimes. Gurarie and Radzihovsky (GR) [45| expanded 
this argument into the vicinity of the BCS-BEC evolu- 
tion with arbitrary winding number k, who finds that 
U£ changes from the Bessel function J? to the modified 
Bessel function I? in the strong coupling regime. Follow- 
ing the procedure in Ref. (45l | , the zero energy solution of 
the BdG equation © with l=(«;+l)/2eZ and E e = 
exhibits 



U( 



(r) = v e (r) = NJi f-^/M 2 -!) e-^, (22a) 



for kn£ > 1 and 



u e (r) = v e (r) =Nh (jyfl-ik^Y 



,-r/i 



(22b) 



for < 1. Here, It(z) is the modified Bessel function 
and TV is the normalization constant. In Eq. (|2"2")l , the 
London approximation is carried out as A+i(r) = Ae 
and A_i(r) = 0. It is obvious that the GR solution in 
Eq. ([25]) coincides with the CdGM solution in Eq. (JT7]) 
in the BCS limit that fc M f;>l. 
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It is seen from Figs. [9(a) and [9(c) that both the 
CdGM and GR solutions within the London approxima- 
tion are in good agreement with the wave functions ob- 
tained from the full numerical calculation, at the BCS 
regime E^/Ep = 4.0 with = 8.6. Here, the de- 
viation around the core can be improved by replacing 
the constant pair potential to A+i(r) = Atanh(r/£ c ), 
where £ c =£/2 is approximately obtained by fitting with 
the self-consistent pair potential, c.f., see the inset of 
Fig. 12(a). In the vicinity of the BCS-BEC topologi- 
cal transition, e.g., E^/Ep = —1.0 with k^ = 0.23 in 
Figs- Sib) andEJd), the GR solution in Eq. flU) repro- 
duces the numerical result of |w^(r)| in contrast to the 
CdGM solution. In this regime, the spatial shape of the 
wave function changes from Jg to the modified Bcsscl 
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FIG. 9: (Color online) (a,b) Wave functions |«f(r)| with £ — 
obtained from the full numerical calculation (solid curves) in 
the k+i = — 1 vortex: (a) Eb/Ep — A.O with fc M £ = 8.6 and (a) 
Eb/Ep — — 1.0 with fc M £ = 0.23. These are compared with the 
GR (dashed lines) and CdGM (dotted curves) solutions with 
the London approximation A(r) = A. The GR and CdGM 
curves are obtained from Eqs. (|22|) and (|17|l . respectively. In 
addition, we plot the CdGM solution obtained from Eq. (|17p 
(dashed-dotted curves), where A(r) = A tanh(r/£ core ) with 
£corc = 2fcp 1 is assumed. (c,d) Wave functions |u^(r)| with £ = 
2 obtained from the full numerical calculation (solid curves) 
in the k+i = +3 vortex: (c) Eb/E-p = 4.0 with k^ = 8.6 and 
(d) E h /E F = -1.0 with fc M £ = 0.23. 



V. CONCLUDING REMARKS 

Here, we have investigated the vortex structures in 
BCS-BEC evolution of p-wave resonant Fermi gases. By 
using the fully microscopic theory based on the BdG 
equation, we have revealed in Sec. Ill the direct relation 
between macroscopic vortex structures and low-energy 
quasiparticle spectra in various vortex states. The low- 
energy spectrum for a single vortex state with vorticity k 
is found to consist of the k branches of the CdGM states 
in addition to the chiral edge state. In particular, it is 
found that in the strong coupling regime, the quantum 
depletion appears in the particle density around the vor- 
tex core, which reflects the shift of the lowest CdGM state 
against the effective pairing interaction. Hence, it is pro- 
posed that the absorption images around vortex cores can 
provide the information on the low-lying CdGM state. 

In addition, in Sec. Ill, we have observed in the k x +ik y 
pairing state with vorticity k that the net angular mo- 
mentum is well estimated as (L z )/H~ (k + 1)N/2 in both 
the BCS and BEC regimes. This turns out to be distinct 
from that in the s-wave case, which can be estimated as 
(L z )/hKi kN/2. Since this deviation originates from the 
orbital angular momentum of the chiral k x + ik y pairing, 
the detection of this deviation can be the direct evidence 
for p-w&ve superfluidity in experiments [52[ . 

While we assume in this paper a rigid boundary condi- 
tion at a certain radius without a harmonic trap poten- 
tial, the Fermi atoms in a realistic system are confined 
by a three-dimensional harmonic trap which may alter 
the conclusions that we present here. However, it is in- 
ferred that such a harmonic trap does not alter the main 
outcome obtained here, as long as the single vortex state 
is assumed. For the CdGM states, this is because the 
typical size £ sa fc~ 1 of the vortex core is much smaller 
than the size of the cloud, of which the Thomas-Fermi 
radius for non-interacting systems is estimated through 
i?TF = \j2E-p /Mlo 2 approximately as an order of 10 /jm, 
where uj is the trap frequency. Hence, the curvature due 
to the trap potential may be much smaller than the spa- 
tial variation of the order parameter around the vortex 
and is ineffective to the results on the CdGM states and 
the depletion of the particle density. It has also been 
demonstrated in our previous work [52j that the trap po- 
tential does not alter the low energy dispersion of the 
edge states in the vortex-free state as well as the net 
angular momentum (L z )/h « N/2. However, such a po- 
tential might work against the stability of the Majorana 
zero modes when the cloud contains many vortices under 
rapid rotation. 

Furthermore, it has been discussed in Sec. IV that the 
zero energy state appears in the BCS phase of an odd 
k vortex state. It is revealed that the lowest eigenener- 
gies exponentially lift from zero to finite values as the 
strong coupling regime is approached. This is because 
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the length scale of the wave function in the CdGM and 
edge state with Eg = becomes comparable to the macro- 
scopic system size and the CdGM and edge states having 
a degenerate energy are hybridized with each other. In 
the extensive range from the weak coupling BCS regime 
to the BCS-BEC topological phase transition point, the 
wave function of the zero energy state is compared with 
the analytic expressions, such as the CdGM solution in 
Eq. dT7J> and the GR solution in Ref. [H|. Here, it is 
demonstrated that the spatial profile turns out to ex- 
hibit not Ji but the modified Bessel function If in the 
strong coupling regime with fc M £ < 1 [45| . 

Finally, it is worth mentioning that the zero energy 
states are unchanged by a slow rotation of the system. 
This can be demonstrated that for axisymmetric vortices 
with arbitrary winding numbers (re+i, kq, re_i) = (re, k + 
1, re + 2), the eigenencrgy Ee(£V) obtained from the BdG 
equation ((5|) under rotation with a frequency f2 yields a 
linear shift from the value without rotation Et(£l = Q) as 



E e (n) = E e (n=o) - £ 



i 



n, 



(23) 



under slow rotation within fl^Ep. This implies that the 
zero energy state with 1= (re + l)/2 € Z remains as it is, 
regardless of S7. In this argument, we neglect the possibil- 
ity on the change of the pair potentials A±i(r), such as 
the penetration of the more vortices. The self-consistent 
pairing field and the existence of the zero energy state in 
a rapidly rotating regime is open for future study. 
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Appendix A: Self-consistent equations 

Here, we describe the details on the derivation of the 
BdG equation §5§ and gap equation First, in order 
to diagonalize the mean- field Hamiltonian in Eq. (TTJ), we 
introduce the unitary transformation to the quasiparticle 
basis with r} u = [rj v , f?J] T , 



(Al) 



where u v is a 2x2 matrix and the matrix elements describe 
the quasiparticle wave functions. They are required to 
satisfy the orthonormal condition, 

J ul{ri)u v <{rx)drx=8 v y, (A2) 



and completeness conditions, ^2 u v (ri)u\,(r2) = 8(r\ — 
r?). Also, the fermion operators r\ v and r\ J obey the 
anti-commutation relation, {^,?7^} = S u y, {rii,,^^} = 
{vtivt'} = 0. The mean-field Hamiltonian in Eq. (fTJ is 
now diagonalized in terms of this basis with the quasi- 
particle energy E v as % = Eq + \ E u ijIt 3 t] u . This di- 
agonalization leads to the Bogoliubov-de Gennes (BdG) 
equation, J dr 2 /S(ri, r 2 )u„(r 2 ) = E v u u (n). Here, it 
is found that the matrix elements of u yield (11)22 = 
(u)*i and (u)i2 = (u)2i, because of the symmetry, 
— fi/C*(ri,r 2 )Ti =/C(ri, r 2 ), where 71,2,3 denote the Pauli 
matrices. Hence, the BdG equation reduces to the equa- 
tion for the quasiparticle wave functions u v = (u)\\ and 
V v = {u)2\, 



dr 2 fc(ri,r 2 ) 



u v (r 2 ) 
v v (r 2 ) 



= E V 



u v (ri) 
v v (ri) 



(A3) 



We then expand the pair potential A(ri, r 2 ) in Eq. ^ 
to the Fourier series with respect to the relative coordi- 
nate ri2=ri — r 2 as, 



/rlk 
(2^ eifc ' ri2A(r ' fc) ' 



(A4) 



with the relative wave vector fc.The Fourier coefficient 
A(r, k) is assumed to be expanded in terms of a p-wave 
channel, 



A(r,k) = T(k) J2 A ™W^« 



(A5) 



m=0,±l 



where we set T(k) = k/k and k±i = ^i(k x ± ik y ) and 
fco = — ik z . We will later see that the pair potential A m (r) 
in the m orbital channel can be expressed in terms of the 
Bogoliubov quasiparticle with the wave function (u u , v u ) 
and the energy E v . 

By substituting Eq. (|A5I) into Eq. (|A4[) and following 
the procedure in Ref. [60( , the off-diagonal element in the 
BdG equation (|A3|) is rewritten as 



(A6) 



A (ri,r 2 ) = ^J2 ^m(r)V^S( ri r 2 ). 



Hence, the BdG equation for the quasiparticlcs under 
pair potentials in m orbital channels are now obtained 
by 



H (r) II(r) \u u (r) 
-II* (r) -H*(r)\ [v v (r) 

m=0,±l - 



= E V 



u v {r) 
v v (r) 



-V m A m (r) 



(A7a) 



, (A7b) 



where the p-wave operator V m is defined as V±\ = 
=F (d x ± id y ) and Vo = d z . The wave functions must sat- 
isfy the orthonormal condition in Eq. (|A2|) 



(A8) 
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Then, in order to derive the gap equation for A m (r), 
we express the Fourier coefficient A(r, k) in Eq. (|A4[) 
with Eqs. (J4J and (|A"Tj) as 



A(r,fc) = y dri2e- iferi2 F(f)$(ri,r 2 ) 



(A9a) 



Eq. (fTB) . from the BdG equation ([5]). Without loss of 
generality, let us consider the situation that k x +ik y pair- 
ing state is majority and k z component is negligible, i.e., 



A +1 (r) = A +1 (r)e i 



(Bla) 



$(ri,r 2 ) = <{r 2 )u„{ ri )f{E v ). (A9b) 

V 

Here, we assume the following p-wave symmetric inter- 
particle interaction, 



V k = ( dr 12 e- ik - r ™V(r) = £ 



m=0,±l 



r(fc)fc m (Aio) 



where g m < is the coupling constant for the m or- 
bital channel of the scattering. Assuming the limit of 
the ri2— >0, the Taylor expansion of the pair function $ 
gives 



$(ri,r 2 ) w $(r,r) 

(Vi- Va)*(n,ra) 



Then, by substituting the interaction V k in Eq. (|A10[) and 
<J> in Eq. (| Al 1|) into the gap equation (|A9|) . one finally 
finds 



A(r,k) 
, .90 



^ 5ro |r(fc)F lim (fe) 



<f>(r, r) 



fc 



r(fc) 



7 (p(.\)_pP0)$( rijPa ) 



fc+i 



(A12) 



Here, the first term should be zero due to the p-wave 
interaction. With the expression in Eq. (|A5j) . one reads 
the gap equations in local form as 



A^r) = ^ fag 



V^-VP )*(ri,r 2 ) 



,(Al3a) 



,(A13b) 



The sum in Eq. (|A13|) denotes the summation for all 
the eigenstates with the positive and negative eigenval- 



A_i(r) = A-^rJe*^ 9 , 



(Bib) 



and Ao(r") = 0. The CdGM solution has been obtained by 
Kopnin and Salomaa 3^ for the negative vortex state, 
where the vorticity is anti-parallel to the chirality of the 
pairing (k = — 1), and by Stone and Chung [47[ for the 
more general case. Also, Volovik [35[ has analytically 
solved the BdG equation ([5]), based on the quasiclassi- 
cal approximation with a quantization rule. As long as 
the zero energy states, the alternative expression was de- 
rived in Refs. 4ol 72l . Here, we expand the expression 



derived in Refs. [3411471 into the more general form which 
is applicable to vortex systems with an arbitrary wind- 
ing number n and the induced component A_i. The 
important consequence obtained here is that the quali- 
tative results on the CdGM state are unchanged by the 
minority A_i component. 

We start with the BdG equation in the cylindri- 
cal coordinate described in Sec. II B. Assuming q = 

& M v/l - sin 2 (a)<fc Al EEE ^2M|/i| 2 , the BdG equation ((SJ) 
reduces to 



M 



dr 

2 



( K + !)(£- £±1) . 



1 dD_ D_ 
k ~27~ 



2 dr 
-T3 



f 2 + 2ME„t 3 



K + l 

2 



n 



(B2) 



where we set D± = D±{r) = A + i(r) ± A_i(r), u v = 
u v (r) = [u„(r), v v (r)] T : and v = (n,l,q). Also, we 
introd uce C m = £2 + \ j r r ~ + k l^ 2 {a) with 
m = \Ji 2 - {k + 1)£+ i*±l)E and 2x2 unit matrix f = 
diag(l,l). Throughout this Appendix, we consider the 
BCS regime with /i>0. 

Following the procedure proposed by Caroli et al. [H| , 
we introduce a radius r c that A(r)=0 for r<r c . Then, 
the BdG equation (|B2[) can be analytically solved if the 
following conditions are assumed: (i) \£\ <C r c fcp <C 
(ii) fc F £>l, (iii) |^|<|M| 2 sin 2 (a), and (iv) |A + i(r)|> 
|A_i(r)|. These conditions restrict to the low-energy 
states in weak coupling BCS regime. 

The solution in Eq. (|B2[) is obtained in the range of 
r < r r as 



Appendix B: Core-bound states in a chiral p-wave 
superfluid 

In this Appendix, we describe the details on the deriva- 
tion of the analytic expression of the CdGM states in 



u v {r) = N 



Je{k + (a)r) 
Je-K-i(k-(a)r) 



where Af is the arbitrary constant and we set 



k± (a) = ku sin(a) ± — ^-r . 



(B3) 



(B4) 
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with Vfj, (a) = sin(a) /M. 

For r > r c , the wave functions arc assumed to consist of 



the Hankel function iJ,„ and the slow functions varying 

(B5) 



over the order of £, <fi(r), as 



Then, Eq. (|B2[) within the conditions (i)-(iii) reduces to 



Tq — + T-2- 
dr 



-IT 3 



MD^(r) 



E v 



(pi(r) = 



n 



D + (r) 



v^a) 2Mv fl (a)r 2 



fi(r), (B6) 



and ip2( r ) T 3 ipl(r). Under the conditions (i)-(iii) de- 
scribed above, the right hand side of Eq. (|B6|) can be 
regarded as the small perturbation. Hence, the regular 
solution within the first order on ij)\,2 is given by 



i(r) = ¥>i 0) (r) +iB ie - x( - r) 



iip 2 (r) 



(B7) 



where ipf^ (r) is the regular solution when the right 
hand side of Eq. (|B6|) is neglected, i.e., <Pi(r) = 
Bie- x ^[l,i] T , and 

M r 

X (r) = — / [A +1 (r')-A_ 1 (r')]dr / . (B8) 



Since |^i j 2(r)|<C!l, Eq. (|B7[) can be also expressed as 

0*01 M 



i(r) ~ Bie-xW 



iijj 2 (r) 



(B9) 



Within Eq. (|B7j) . one can find the solution of Eq. (|B6 
as t/'i( r ') = — ip2( r ) = ip(r), where 



u M (a) 2Mv^(a)r' 



D + (r>) e-*p 
v-p Mvfj,(a)r' 



In order to obtain the solution of the BdG equation 
(|B2[) . the wave functions in Eq. (|B3|) for r < r c and 
Eq. ()B5[) for r > r c , are now matched at r = r c . Be- 
cause of the condition (i), \£\ <C r c /cp, making use of the 

asymptotic forms of Je(z) and i/^ : (z) in z> \£\, the 
wave functions for r < r c in Eq. (|B3[) is rewritten as 



u„ w TV- 



2M 



£ 2 - \ 21+1 
cos K+r H 2- 7T 

4 2£ v + 1 



cos fc 



2fc_r 



(Bll) 



with i> M = Vfj,(a) and fc± = fc±(a). Also, Eq. (|F35|) with 
Eq. flBSJ} for rfc F >r c fc F >|^| is 



with 



/ 2M 



,-xM 



, (B12) 



?7±(r) = A/w^r + 



m 2 - | 2m + 1 
2Mv~r ~ 4 



7r±V(r). (B13) 



Then, the coefficients i?!^ are determined so as to 
match two expressions of ttjy(r) in Eqs. (jBTT|) and (|B12|) 

at r = r r as 



(B14) 



By comparing with Eqs. (|B1 1|) and (|B12|) . one can obtain 
the expression of ip as 



u M (a) 2M^(a)r 



(B15) 



In the same way, the another expression is obtained from 
Eqs. (|B"TT|) and ([Bl"2j) 



(«+!)(£_ «±1) 

-r - ' 



u M (a) 2Mw )1 (a)r 
-^(m-^ + K + l)+7-mr, (B16) 



where n£Z. The expressions on ^>(r) in Eqs. (|B15j) and 
(|B1G|) becomes identical when 7 satisfies 



7 = g 1 m 



K + 1 \ 7T 

— + 2 n - 



(B17) 



The alternative expressions of ^f(r) in Eq. (jBlOl) and 
Eq. (|B15j) with Eq. ()B17|) should be identical at r = r c . 
Hence, we finally obtain the eigenvalue of the BdG equa- 
tion (|A"7l) , 



E v = 
where 

oj = 



k + i\ ( k + r. 

Wo + [ n — ) sin(a)wi, (B18) 



«+iA+i(r') - n-i A^ir*) c _ 3x(r , )dr , 



kj?T f 



nk n 



2M / e- 2x(r ' ] dr' 
Jo 



(B19a) 



(B19b) 
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The resulting eigenvalue in Eq. (|B18|) reproduces the re- 
sults derived in Res. [U, H3] and the qualitative proper- 
ties, such as the appearance of the ZES, turn out to be 
unchanged by the minority component A_i. 

The eigenfunction is then obtained from Eqs. (|B12|) - 
(|BT7|) with k± ( a) in Eq. (|B4j) as 



= TV 



Je(k + (a)r) 
Ji- K -i(k-(a)r) 



cxp f [A +1 (r') - A_!(r')] dr'j .(B20) 



To estimate the order of the energy scale of wo,i, let us 
consider the simplest case of A(r), that is, A + i(r) = 
Atanh(r/£) and A_i(r)=0. In this situation, one finds 
wi « ^A and w ~ with J°° e _2x ^ r W = £. Hence, 
the eigenvalue E v consists of two different energy scales, 
such as A and A 2 /Ep. The expression in the case of 
a singly quantized vortex (k = ±1) coincides to that in 
Refs. [HI and [13], which yields the zero energy modes 
with £ = for k = — 1 and l — — \ for k = +1. 
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